#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int main() 
{
/* -------------------------------------------------------
   -------- Reads the NPHASE matrix in CRS format --------
   -------- Created by Igor A. Bolotnov, 03/01/2009 ------
   -------------------------------------------------------
*/
    FILE *fpin, *fpout;
    int nrow, nvar; // number of rows (nrow) in the block matrix and number of unknown variables (nvar) within each submatrix 
    int irow, icol, icolen, icolst, n, m, irow2;
    int blksize, blksize2, xindex, nbrow;
    double product;
    double *a_matrix, *b_vector, *x_vector, *error_vector; // thos are matrix A; r.h.s. vector B and solution vector X
    int *row_ptr, *col_ind; 	

    printf("This code reads an NPHASE matrix in CRS format and verifies that the solution is correct by computing Ax = B \n");

    // read case data
    fpin = fopen("case.dat0", "r");
    fscanf(fpin, "%d ", &nrow);
    fscanf(fpin, "%d ", &nvar);
    close(fpin);

    printf(" nrow = %d, nvar = %d \n", nrow, nvar) ;

    blksize = nvar;
    blksize2 = nvar*nvar;

    row_ptr = calloc(nrow+1, sizeof(int)) ;

    // read row_ptr data
    fpin=fopen("row_ptr.dat","r");
    for(irow=0; irow<=nrow; ++irow)
    {
	fscanf(fpin,"%d",(row_ptr+irow));
	printf("irow=%d row_ptr=%d \n",irow,*(row_ptr+irow));
    }
    close(fpin);

    // read column index data
    col_ind = calloc(*(row_ptr+nrow)+1, sizeof(int)) ;
    fpin=fopen("col_ind.dat","r");
    icolst= *(row_ptr+0) ;
    icolen= *(row_ptr+nrow-1+1)-1;
    for(icol=icolst; icol<=icolen; ++icol)
    {
	fscanf(fpin,"%d",(col_ind+icol));
        // printf("%d \n",*(col_ind+icol));
    }
    close(fpin);

    printf("col_ind.dat file has been read.  \n");

    // read matrix A data
    a_matrix = calloc(blksize2* *(row_ptr+nrow-1+1), sizeof(double));

    printf("The matrix size is %d \n", blksize2* *(row_ptr+nrow-1+1)) ;

    fpin=fopen("matrix_coef.dat","r") ;
    for(irow=0; irow<=nrow-1; ++irow)
    {
	icolst= *(row_ptr+irow) ;
	icolen= *(row_ptr+irow+1)-1;
	for(icol=icolst; icol<=icolen; ++icol)
	{
	    printf(" irow=%d col_ind=%d \n",irow,*(col_ind+icol));
	    for(m=0; m<=blksize-1; ++m) 
	    {
		for(n=0; n<=blksize-1; ++n) 
		{
		    fscanf(fpin,"%le",(a_matrix+icol*blksize2+n+m*blksize));
		    printf(" %1e ",*(a_matrix+icol*blksize2+n+m*blksize));
		}
		printf("\n");
	    }
	}
    }
    close(fpin);

    printf("The matrix A has been read. \n") ;

    printf("sizeof: %d", sizeof(a_matrix));

    // read RHS Vector B
    b_vector = calloc(blksize* nrow, sizeof(double)) ;

    fpin=fopen("rhs_value.dat","r");
    for(irow=0; irow<=nrow-1; ++irow)
    {
	for(m=0; m<=blksize-1; ++m) 
	{
	    fscanf(fpin,"%le",(b_vector+irow*blksize+m));
	}
    }
    close(fpin);

    printf("The vector B has been read. \n") ;

    // read X Solution Vector
    x_vector = calloc(blksize* nrow, sizeof(double)) ;

    fpin=fopen("x_value.dat","r");
    for(irow=0; irow<=nrow-1; ++irow)
    {
	printf(" irow=%d \n",irow);
	for(m=0; m<=blksize-1; ++m) 
	{
	    fscanf(fpin,"%le",(x_vector+irow*blksize+m));
	    printf(" %1e ",*(x_vector+irow*blksize+m));
	}
	printf("\n");
    }
    close(fpin);

    printf("The NPHASE linear solver solution vector X has been read. \n") ;

    // Now we have to check that solution satisfies the equation, this way we we certain that we properly have read
    // and understood the matrix structure

    // Vector error_vector is A*X-B

    error_vector = calloc(blksize* nrow, sizeof(double)) ;

    //	printf("b_vector[0] = %le \n", *(b_vector+0)) ;

    // Loop over rows in the block matrix:
    for (irow=0; irow<=nrow-1; ++irow) 
    {
	// loop over rows in the submatrices:
	for (m=0; m<=blksize-1; ++m) 
	{
            // here we have to accumulate the error_vector value:
	    nbrow = irow*blksize + m ;  // number of the vector component being computed (matrix row)
	    *(error_vector+nbrow) = (double) -*(b_vector+nbrow); 		// first step: error_vector = -B
	    printf("Processing component N%d \n", nbrow) ;
            // loop over blocks in the row "irow"
	    icolst= *(row_ptr+irow) ;			// start column location in A
	    icolen= *(row_ptr+irow+1)-1;			// end column location in A
	    for(icol=icolst; icol<=icolen; ++icol)
	    { 	// go over the block columns (in the block row "irow")
		// the column number in the A matrix which corresponds to the row number in vector X. We use this value as a running variable as we multiply the row nbrow of A by vector X  
		xindex = *(col_ind+icol)*blksize ;  	
		for(n=0; n<=blksize-1; ++n) 
		{   // go over the submatrix columns
		    // Note that X row has to be equal to the A column (since we are multiplying a vector by a matrix !)
		    product = (double) *(a_matrix+icol*blksize2+n+m*blksize) * *(x_vector+xindex+n) ;
		    *(error_vector+nbrow) += product ;	// here error_vector = B - AX
		    if (nbrow == 5) 
		    {
			printf("Accumulating the product of the matrix component ((icol = %d)column %d, row %d) and vector X component %d in the error_vector component %d \n",
			       icol, *(col_ind+icol)*blksize+n, nbrow, xindex+n, nbrow);
			
			printf("icol = %d, n = %d, col_ind = %d, Bindex = %d, a_matrix = %le, b_vector = %le, x_vector = %le error_vector = %le \n", icol, 
			       n,*(col_ind+icol), xindex,  *(a_matrix+icol*blksize2+n+m*blksize), *(b_vector+xindex+n),  *(x_vector+xindex+n), *(error_vector+nbrow)) ;
		    } 	//if
		} 	// n
	    } 	// icol
	}  	// for m
    } 	// for irow
    
    // Print out the error_vector vector:
    fpout=fopen("zero_check.dat","w");
    for(irow=0; irow<=nrow-1; ++irow)
    {
	for(m=0; m<=blksize-1; ++m) 
	{
	    fprintf(fpout,"%le\n",*(error_vector+irow*blksize+m));
	}
    }
    close(fpout);

    printf("     ... done \n");

    free(row_ptr);
    free(col_ind);
    free(a_matrix);
    free(b_vector); 
    free(x_vector);
    free(error_vector);

    return 0 ;
}
